function [averageBetas maskedBetas averageBetasMatrix stdBetasMatrix] = extractROImask(subjID, f_thresh, FFAcontrast, analysesDir)
%AdaptID_v1_2_Faces: Call the AdaptID_v1_2  driver for face stimuli
%           [params data stimSet] = AdaptID_v1_2_Faces(expInfo)
%
%           AdaptID_v1_2_Faces(expInfo) initializes the driver for the
%   AdaptID_v1_2 experiment. This file should not be called directly.
%   Instead, it should be called from the AdaptID_v1_2() function.
%           The driver reads the experimental parameters, sets up the
%   screen, load the stimuli, and run the experiment. It also collects and
%   saves the relevant data.
%
%
%   Inputs
%   -------
%   expInfo: contains all the experimental information necessary for
%           execution, including the experimental protocol to be ran
%
%   Outputs
%   -------
%   params: Experimental parameters, read from an external file specified
%           on the protocol file
%   data: Data collected during execution, including subject responses,
%           reaction times, and all presentation times
%   stimSet: Information about the stimuli that were presented
%
%   _______________________________________________
%   by Marcelo G Mattar (09/04/2012)
%   mattar@sas.upenn.edu


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% CHECK INPUTS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin < 4
    analysesDir = '/Volumes/cluster/jet/mattar/AdaptID/Analyses/';
end

if nargin < 3
    FFAcontrast = 2; %1 for Faces-Objects, 2 for Faces-Scenes
end

if nargin < 2
    f_thresh = -Inf;
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% BASIC DEFINITIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% Define some basic directories and files
subjectDir = [analysesDir subjID '/'];
FOBSdir = [subjectDir 'FOBSLoc.feat/'];


% Retrieve the analysis directories
FaceDirectories = struct2cell(dir([subjectDir 'Face*.feat']));
FaceDirectories = FaceDirectories(1,:)';
numFaceRuns = length(FaceDirectories);


averageBetas = cell(numFaceRuns+2,8);
averageBetas{1,1} = subjID;

averageBetas{2,1} = 'Run #1';
averageBetas{3,1} = 'Run #2';
averageBetas{4,1} = 'Run #3';
if numFaceRuns == 4
    averageBetas{5,1} = 'Run #4';
end
averageBetas{numFaceRuns+2,1} = 'Average';

averageBetas{1,2} = '1-back';
averageBetas{1,3} = '2-back';
averageBetas{1,4} = '3-back';
averageBetas{1,5} = '4-back';
averageBetas{1,6} = '5-back';
averageBetas{1,7} = '6-back';
averageBetas{1,8} = 'F-test';

averageBetasMatrix = zeros(numFaceRuns,6);
stdBetasMatrix = zeros(numFaceRuns,6);


% Pre-allocate some important variables
feat_cope = cell(numFaceRuns,6);
feat_varcope = cell(numFaceRuns,6);
feat_tstat = cell(numFaceRuns,6);

gfeat_cope = cell(1,6);
gfeat_varcope = cell(1,6);
gfeat_tstat = cell(1,6);

leaveOneOut_cope = cell(numFaceRuns,1);
leaveOneOut_varcope = cell(numFaceRuns,1);
leaveOneOut_tstat = cell(numFaceRuns,1);

maskedBetas = cell(numFaceRuns,6);



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% LOAD FFA MASK
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if exist([FOBSdir 'stats/lFFA' num2str(FFAcontrast) 'mask.nii.gz'],'file')
    lFFAmask = loadniigz([FOBSdir 'stats/lFFA' num2str(FFAcontrast) 'mask.nii.gz']);
else
    possibleContrasts = [1 2];
    lFFAmask = loadniigz([FOBSdir 'stats/lFFA' num2str(possibleContrasts(~(FFAcontrast == [1 2]))) 'mask.nii.gz']);
end

if exist([FOBSdir 'stats/lFFA' num2str(FFAcontrast) 'mask.nii.gz'],'file')
    rFFAmask = loadniigz([FOBSdir 'stats/rFFA' num2str(FFAcontrast) 'mask.nii.gz']);
else
    possibleContrasts = [1 2];
    rFFAmask = loadniigz([FOBSdir 'stats/rFFA' num2str(possibleContrasts(~(FFAcontrast == [1 2]))) 'mask.nii.gz']);
end

FFAmask = lFFAmask + rFFAmask;
FFAmask = logical(FFAmask);






%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% LOAD ADAPT DATA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Extract Adaptation copes, varcopes and tstats for each of the Face runs
for runIndx = 1:numFaceRuns
    for adaptIndx = 1:6
        feat_cope{runIndx,adaptIndx} = loadniigz([subjectDir FaceDirectories{runIndx} '/reg_standard/stats/cope' num2str(adaptIndx+1) '.nii.gz']);
        %feat_varcope{runIndx,adaptIndx} = loadniigz([subjectDir FaceDirectories{runIndx} '/reg_standard/stats/varcope' num2str(adaptIndx+1) '.nii.gz']);
        %feat_tstat{runIndx,adaptIndx} = feat_cope{runIndx,adaptIndx} ./ sqrt(feat_varcope{runIndx,adaptIndx});
    end
end

FstatsMatrix = loadniigz([subjectDir 'Face_Ftest.gfeat/cope1.feat/stats/fstat1.nii.gz']);
thisMask = and((FstatsMatrix > f_thresh), FFAmask);

for runIndx = 1:numFaceRuns
    for adaptIndx = 1:6
        maskedBetas{runIndx,adaptIndx} = feat_cope{runIndx,adaptIndx}(thisMask(:));
        averageBetasMatrix(runIndx,adaptIndx) = mean(maskedBetas{runIndx,adaptIndx});
        stdBetasMatrix(runIndx,adaptIndx) = std(maskedBetas{runIndx,adaptIndx});
    end
end


